Introduction à l’analyse des données

1 Markdown 101

Markdown est un langage de balisage et peut être utilisé pour donner aux documents texte tels que R Markdown une structure comme word ou latex. La syntaxe est assez simple.

quelques exemples (tiré de https://www.markdownguide.org/basic-syntax):

# Heading level 1   
Titre de section


## Heading level 2  
Titre de sous-section

### Heading level 3
Titre de sous-sous section

etc..

__un texte en gras__

_un texte en italique_

~~un text barré~~

ie

un texte en gras

un texte en italique

un text barré

Pour plus d’informations, consultez par exemple le lien de rstudio.

1.1 LaTeX

R Markdown est encore plus puissant, on peut utilisier \(\LaTeX\), soit inline comme \(y_i=\beta_0 + x_i\beta_1\) ou comme equation: \[ \hat{\beta}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{y} \] vous pouvez même inclure des fichiers .tex entiers.

1.2 Figures

Imuages (ce n’est pas un plot) peuvent également être facilement incorporés avec ![](la_meilleure_image.png)

1.3 Knitting

Pour les documents de base, le yaml est ce qui est important:

rmdformats::robobook c’est seulement pour rendre le ficher légèrement plus joli

1.4 Pourquoi?

  • Projets pour votre thèse de baccalauréat (Peut être entièrement reproductible)
  • après votre baccalauréat
    • Industrie-1: Vous ferez beaucoup de rapports, souvent répétitifs, ce qui n’est pas très fun.. mais vous pouvez automatiser la plupart avec markdown
    • Industrie-2: Particulierement en tant que junior, vous aurez besoin de faire des recherches, Rmarkdown est parfait pour les présenter
    • Université: Si vous faites une maitrise - Vous devrez probablement coder et remettre votre code
    • Vous pouvez également utiliser python à partir des chunks de rmarkdown

2 Manipulation des Données

Si vous travaillez avec des données, vous passerez beaucoup de temps à préparer vos données avant de modeliser (ensembles de données “propres” n’existent généralement pas dans le monde réel). Dans le cadre de ce cours, nous n’examinerons que les principes de base des “tidy data”.

L’idée générale est de mettre de l’ordre dans nos données de façon à ce que nous ayons :

  • une observation par ligne
  • Une information par cellule
  • Un type d’information par colonne

Pour plus il y a le résumé du papier de Hadley Wickham, vous pouvez consulter ce lien

Dans ce qui suit, nous envisagerons toujours plusieurs façons de résoudre un problème. Vous êtes libre de choisir celle que vous correspond le plus, bien que je recommande personnellement l’écosystème “tidyverse”.

Pour des raisons esthétiques, il est souvent recommandé d’utiliser l’opérateur pipe (%>%). Sous R-Studio avec Cmd-Shift-M ou Ctrl-Shift-M. Ce qu’il fait fondamentalement est de prendre la ligne précédente comme premier argument dans la ligne courante.

# meme chose
my_numbers <- c(0,1,2,3)

print(sum(my_numbers))
## [1] 6
print(my_numbers %>% sum())
## [1] 6

cette logique simplifiera grandement vos flux de travail et rendra le code plus lisible.

2.1 SQL-Like

SQL (Structured Query Language) - votre meilleur ami si vous travaillez avec des données dans une entreprise, mais également utile à plus petite échelle (aka. pour vos projets).

L’exercice est basé sur l’article : “Where is the Land of Opportunity ? The Geography of Intergenerational Mobility in the United States” (Raj Chetty, Nathaniel Hendren, Patrick Kline et Emmanuel Saez ; Quarterly Journal of Economics 129(4) : 1553-1623, 2014).Il fait partie d’un plus grand projet disponible à l’adresse suivante : http://www.equality-of-opportunity.org/.

Nous commençons par charger quelques données échantillons. Notez que “kable” dans la dernière ligne rend seulement les graphiques jolis.

library(haven) # si le format est .dta
library(readxl) # se le format est .xls / .xlsx

mobility_data <- read_dta("data/mobility_across_cz.dta")
codebook <- read_excel("data/codebook.xls")

mobility_data %>% 
  head() %>% 
  kable()
index cz czname e_rank_b cs_race_bla cs_race_theil_2000 cs00_seg_inc cs00_seg_inc_pov25 cs00_seg_inc_aff75 frac_traveltime_lt15 hhinc00 gini inc_share_1perc gini99 frac_middleclass taxrate subcty_total_expenditure_pc tax_st_diff_top20 eitc_exposure ccd_exp_tot ccd_pup_tch_ratio score_r dropout_r num_inst_pc tuition gradrate_r cs_labforce cs_elf_ind_man d_tradeusch_pw_1990 frac_worked1416 mig_inflow mig_outflow cs_born_foreign scap_ski90pcm rel_tot crime_violent cs_fam_wkidsinglemom cs_divorced cs_married
0 100 Johnson City 43.58841 0.0693926 0.0576880 0.0408645 -0.0171138 0.0266398 0.3844612 38209.54 0.4193436 5.695437 0.2528637 0.5840611 0.0158913 1404.822 1.0794937 0.0000780 6.577284 17.96812 NA NA 0.0554321 1032.0113 -0.2210098 0.6529477 0.0600338 NA 0.0148774 0.0417411 0.0559070 0.1604878 NA 0.5358987 0.0482720 0.2421846 0.1392865 0.5491035
1 200 Morristown 39.64703 0.4453178 0.0749721 -0.0258754 0.0352996 0.0089049 0.3543061 25123.34 0.4867183 11.816536 0.3926104 0.3559214 0.0257373 1568.710 0.0198662 -0.0128659 4.405095 17.79943 -16.015344 0.0561887 0.0379247 970.0455 -0.0416078 0.5273348 0.2086723 0.7670213 0.0223866 0.0181590 0.0137850 0.0232165 -0.6127989 0.6469429 0.0316044 0.3499764 0.1039057 0.5426201
2 301 Middlesborough 46.93892 0.0151223 0.1100580 0.0954185 0.0139200 0.0527287 0.4793166 32701.43 0.3071970 11.810008 0.2077348 0.6675747 0.0116976 3115.946 -0.0046168 0.0343010 6.609784 17.85102 -9.137633 0.0246227 0.0634999 1395.0190 -0.1716805 0.7066224 0.0362908 0.0274560 0.0274638 0.0309692 0.0421808 0.1318125 -2.0635338 0.4222386 0.0511852 0.1831544 0.1402021 0.6224286
3 302 Knoxville NA -0.0004112 -0.0044561 0.0104628 0.0545056 0.0076278 0.6607098 36523.62 0.2411462 NA NA NA 0.0216987 2513.664 0.2094634 1.9299725 7.892078 12.60501 NA NA NA NA NA 0.6233307 0.0685092 0.0173301 NA NA NA 0.0252765 2.3432735 0.6877957 0.0371077 0.1290979 0.1127880 0.7102951
4 401 Winston-Salem 43.22648 -0.0193413 0.2634850 0.0288748 0.0478488 0.0685591 0.4064568 30403.07 0.4511340 9.971601 0.3481184 0.5376577 0.0451547 2137.595 0.0028042 -0.0035268 7.054153 20.89541 -10.453489 NA 0.0418891 2523.8348 0.1363400 0.6357111 0.1408433 0.8333290 0.0239467 0.0580380 0.0682463 0.1903379 -0.6936674 0.3871659 0.0499824 0.2215046 0.1146979 0.5695345
5 402 Martinsville 53.63568 0.0056276 0.0069904 0.0220931 0.0410883 0.0172575 0.5862997 29943.32 0.2565523 4.503661 0.1837683 0.7311252 0.0266855 2756.848 0.3987184 -0.0275933 4.544135 18.19494 5.789057 -0.0153886 NA NA NA 0.6121494 0.1984681 0.5390288 0.0377694 0.0170314 0.0643845 0.0236150 0.3655432 0.8407629 0.0321638 0.1507503 0.1158269 0.6865350

Nous voyons qu’il y a des données manquantes et que les variables explicatives sont normalement numériques. Nous nous en occuperons plus tard. Pour une explication des variables, considérez le fichier codebook.xls.

select

En général, un jeu de données comporte de nombreuses colonnes. Mais parfois, nous n’avons besoin de considérer que certaines d’entre elles. La logique “select” vous permet de choisir un sous-ensemble des colonnes de données.

Il existe de nombreuses façons de le faire. Ci-dessous la démonstration pour sélectionner une seule colonne “e_rank_b” et une collection de colonnes c(“e_rank_b”, “taxrate”, “crime_violent”). C’est quoi la difference?

# Une seule colonne
facon_1 <- mobility_data$e_rank_b
facon_2 <- mobility_data[,'e_rank_b']
facon_3 <- mobility_data %>% select(e_rank_b)

# Plusieurs
facon_1_mul <- mobility_data[, c("e_rank_b", "taxrate", "crime_violent")]
facon_2_mul <- mobility_data %>% select(e_rank_b, taxrate, crime_violent)


facon_1 %>% head() 
## [1] 43.58841 39.64703 46.93892       NA 43.22648 53.63568
facon_1_mul %>% head()
## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     39.6  0.0257        0.0316
## 3     46.9  0.0117        0.0512
## 4     NA    0.0217        0.0371
## 5     43.2  0.0452        0.0500
## 6     53.6  0.0267        0.0322
facon_2 %>% head() 
## # A tibble: 6 × 1
##   e_rank_b
##      <dbl>
## 1     43.6
## 2     39.6
## 3     46.9
## 4     NA  
## 5     43.2
## 6     53.6
facon_2_mul %>% head() 
## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     39.6  0.0257        0.0316
## 3     46.9  0.0117        0.0512
## 4     NA    0.0217        0.0371
## 5     43.2  0.0452        0.0500
## 6     53.6  0.0267        0.0322

if-condition (filter)

Les commandes de “select” agissent sur les colonnes, les équivalents pour les lignes sont les “if-conditions” dans R, on appel ca les conditions “filter”. En reprenant l’exemple ci-dessus. Nous ne voulons que les observations où le e_rank_b est supérieur à 40. Quelle différence peut-on constater entre les deux méthodes ?

facon_1 <- mobility_data[mobility_data$e_rank_b > 40,
                         c("e_rank_b", "taxrate", "crime_violent")]

facon_2 <- mobility_data %>% 
  select(e_rank_b, taxrate, crime_violent) %>% 
  filter(e_rank_b > 40)

facon_1 %>% head()
## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     46.9  0.0117        0.0512
## 3     NA   NA            NA     
## 4     43.2  0.0452        0.0500
## 5     53.6  0.0267        0.0322
## 6     45.6  0.0412        0.0267
facon_2 %>% head()
## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     46.9  0.0117        0.0512
## 3     43.2  0.0452        0.0500
## 4     53.6  0.0267        0.0322
## 5     45.6  0.0412        0.0267
## 6     48.8  0.0222        0.0286

groupby

Parfois, nous sommes intéressés par des données agrégées. Avec la commande group_by, nous pouvons facilement créer des statistiques sommaires pour des sous-groupes de données.

Il existe des façons de faire exactement cela dans base R, mais ici nous nous concentrons sur la logique de dplyr, qui est plus intuitive.

# Calculer la moyenne de "e_rank_b" pour chaque ville

group_stat <- mobility_data %>% # Choisir jeu des données
  select(czname, e_rank_b) %>% # Choisir les colonnes
  group_by(czname) %>% # Specifier les groups
  summarise(
    mean_per_city = mean(e_rank_b, na.rm=T) # Specifier la statistique d'aggreagtion
  )

group_stat %>% head()
## # A tibble: 6 × 2
##   czname     mean_per_city
##   <chr>              <dbl>
## 1 Aberdeen            45.6
## 2 Aiken               38.5
## 3 Ainsworth           48.6
## 4 Albany              52.9
## 5 Alexandria          46.6
## 6 Allentown           61.3

(mutate)

Souvent, nous devons transformer les variables pour obtenir des informations significatives de notre analyse. Cela peut être fait facilement avec la logique “mutate”. L’utilisation de condtions if-else pour créer de nouvelles variables est extrêmement utile. Par exemple, nous pouvons diviser la mobilité sociale (e_rank_b) en deux catégories : faible et élevée.

data_tmp <- mobility_data
data_tmp <- data_tmp[!is.na(data_tmp$e_rank_b), ]
data_tmp$binary_mobility <- ifelse(data_tmp$e_rank_b > 45, 'High', 'Low')

ex_mutate_2 <- mobility_data %>% 
  filter(!is.na(e_rank_b)) %>% # PQ? 
  mutate(binary_mobility = if_else(e_rank_b > 45, 'High', 'Low')) %>% 
  select(binary_mobility, e_rank_b)


data_tmp[,c('binary_mobility', 'e_rank_b')] %>% head()
## # A tibble: 6 × 2
##   binary_mobility e_rank_b
##   <chr>              <dbl>
## 1 Low                 43.6
## 2 Low                 39.6
## 3 High                46.9
## 4 Low                 43.2
## 5 High                53.6
## 6 Low                 38.1
ex_mutate_2 %>% head()
## # A tibble: 6 × 2
##   binary_mobility e_rank_b
##   <chr>              <dbl>
## 1 Low                 43.6
## 2 Low                 39.6
## 3 High                46.9
## 4 Low                 43.2
## 5 High                53.6
## 6 Low                 38.1

2.2 Exercises

  1. Transformer la variable “e_rank_b” en “inférieure à 50” et “supérieure à 50”. Calculez la moyenne de “crime_violent” pour chaque groupe. Quelles sont vos conclusions?
  2. Transformer la variable “frac_middleclass” en “inférieure à 0.5” et “supérieure à 0.5”. Pareil pour la variable “eitc_exposure”. Calculer la moyenne du taux d’abandon scolaire au lycée “dropout_r” pour chaque combinaison des deux nouvelles variables. Conseil : utilisez group_by
  3. Proposer une façon de traiter les valeurs “NA”. Pourquoi pensez-vous que c’est une bonne idée ?

3 Analyse Visuelle

Il est toujours bon d’avoir une compréhension intuitive des données en utilisant des graphiques. Ici, nous essaierons de réaliser la plupart des graphiques à la main (il existe des logiciels qui peuvent faire le travail pour vous, mais il est préférable de le comprendre d’abord complètement).

3.1 Plots

Les plots/figures peuvent facilement être intégrées dans markdown. Tout comme une sortie de tableau, on peut les structurer légèrement. Par exemple, nous pouvons ajouter une légende ou ajouter plusieurs graphiques côte à côte.

Par exemple, commençons par un histogramme:

# Histogram simple
hist(mobility_data$e_rank_b)

?hist

?hist
plot(density(data_tmp$e_rank_b)) # Pourquoi changer les données?
plot(density(data_tmp$e_rank_b), col='red')

Il existe aussi un package sympa de l’écosystème tidyverse (ggplot2)

plot_1 <- mobility_data %>% 
  ggplot() + 
  geom_histogram(aes(x=e_rank_b))

plot_2 <- mobility_data %>% 
  ggplot() + 
  geom_histogram(aes(x=crime_violent))

grid.arrange(plot_1, plot_2, ncol=2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 19 rows containing non-finite values (stat_bin).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 17 rows containing non-finite values (stat_bin).

Nous pouvons également construire des figures avec plusieurs informations. Par exemple, un histogramme avec une ligne de densité. (Attention toutefois aux axes!)

# Base R
hist(data_tmp$e_rank_b, probability = T)
lines(density(data_tmp$e_rank_b), col='red')

# GGplot
data_tmp %>% 
  ggplot() + 
  geom_histogram(aes(x=e_rank_b, y=..density..)) + 
  geom_density(aes(x=e_rank_b), color='red')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

3.2 Exercises

Conseil : si vous n’arrivez pas à trouver la bonne commande pour une figure, les solutions, utilisez google!

  1. Montrer la relation entre le taux d’obtention d’un diplôme universitaire et la mobilité socialeà l’aide d’un diagramme de dispersion (scatterplot) et d’une ligne de régression. Conseil: vous pouvez trouver une solution si vous googlez “r plot with regression line”

  2. Existe-t-il une relation entre les dépenses des gouvernements locaux par habitant (subcty_total_expenditure_pc) et le taux d’imposition (taxrate)? Pouvez-vous utiliser les données directement ? Conseil : pour mettre les données à l’échelle, vous pouvez toujours essayer d’utiliser log, sqrt, etc.

  3. Montrez graphiquement qu’il existe une ville où la criminalité violente est exceptionnellement élevée (‘crime_violent’). Que pouvez-vous faire avec de telles données (outlier ou pas)?

4 Analyse descriptive

L’analyse graphique est utile mais se complique avec plus que deux variables. En général, nous évaluons également les hypothèses que nous avons dérivées des graphiques à l’aide d’une sorte de test.

4.1 Tables

Les tableaux peuvent facilement être incorporés dans R. Vous pouvez soit les imprimer directement, soit utiliser kable (du paquet knitr - plus jolie). Bien sûr, nous pouvons également combiner ce que nous avons vu dans les sections précédentes (en particulier la manipulation des données)

# Impression de base
mobility_data %>% head()
## # A tibble: 6 × 39
##   index    cz czname         e_rank_b cs_race_bla cs_race_theil_20… cs00_seg_inc
##   <dbl> <dbl> <chr>             <dbl>       <dbl>             <dbl>        <dbl>
## 1     0   100 Johnson City       43.6    0.0694             0.0577        0.0409
## 2     1   200 Morristown         39.6    0.445              0.0750       -0.0259
## 3     2   301 Middlesborough     46.9    0.0151             0.110         0.0954
## 4     3   302 Knoxville          NA     -0.000411          -0.00446       0.0105
## 5     4   401 Winston-Salem      43.2   -0.0193             0.263         0.0289
## 6     5   402 Martinsville       53.6    0.00563            0.00699       0.0221
## # … with 32 more variables: cs00_seg_inc_pov25 <dbl>, cs00_seg_inc_aff75 <dbl>,
## #   frac_traveltime_lt15 <dbl>, hhinc00 <dbl>, gini <dbl>,
## #   inc_share_1perc <dbl>, gini99 <dbl>, frac_middleclass <dbl>, taxrate <dbl>,
## #   subcty_total_expenditure_pc <dbl>, tax_st_diff_top20 <dbl>,
## #   eitc_exposure <dbl>, ccd_exp_tot <dbl>, ccd_pup_tch_ratio <dbl>,
## #   score_r <dbl>, dropout_r <dbl>, num_inst_pc <dbl>, tuition <dbl>,
## #   gradrate_r <dbl>, cs_labforce <dbl>, cs_elf_ind_man <dbl>, …
# Avec Kable
mobility_data %>% head() %>% kable() # c'est comme kable(head(mobilitiy_data))
index cz czname e_rank_b cs_race_bla cs_race_theil_2000 cs00_seg_inc cs00_seg_inc_pov25 cs00_seg_inc_aff75 frac_traveltime_lt15 hhinc00 gini inc_share_1perc gini99 frac_middleclass taxrate subcty_total_expenditure_pc tax_st_diff_top20 eitc_exposure ccd_exp_tot ccd_pup_tch_ratio score_r dropout_r num_inst_pc tuition gradrate_r cs_labforce cs_elf_ind_man d_tradeusch_pw_1990 frac_worked1416 mig_inflow mig_outflow cs_born_foreign scap_ski90pcm rel_tot crime_violent cs_fam_wkidsinglemom cs_divorced cs_married
0 100 Johnson City 43.58841 0.0693926 0.0576880 0.0408645 -0.0171138 0.0266398 0.3844612 38209.54 0.4193436 5.695437 0.2528637 0.5840611 0.0158913 1404.822 1.0794937 0.0000780 6.577284 17.96812 NA NA 0.0554321 1032.0113 -0.2210098 0.6529477 0.0600338 NA 0.0148774 0.0417411 0.0559070 0.1604878 NA 0.5358987 0.0482720 0.2421846 0.1392865 0.5491035
1 200 Morristown 39.64703 0.4453178 0.0749721 -0.0258754 0.0352996 0.0089049 0.3543061 25123.34 0.4867183 11.816536 0.3926104 0.3559214 0.0257373 1568.710 0.0198662 -0.0128659 4.405095 17.79943 -16.015344 0.0561887 0.0379247 970.0455 -0.0416078 0.5273348 0.2086723 0.7670213 0.0223866 0.0181590 0.0137850 0.0232165 -0.6127989 0.6469429 0.0316044 0.3499764 0.1039057 0.5426201
2 301 Middlesborough 46.93892 0.0151223 0.1100580 0.0954185 0.0139200 0.0527287 0.4793166 32701.43 0.3071970 11.810008 0.2077348 0.6675747 0.0116976 3115.946 -0.0046168 0.0343010 6.609784 17.85102 -9.137633 0.0246227 0.0634999 1395.0190 -0.1716805 0.7066224 0.0362908 0.0274560 0.0274638 0.0309692 0.0421808 0.1318125 -2.0635338 0.4222386 0.0511852 0.1831544 0.1402021 0.6224286
3 302 Knoxville NA -0.0004112 -0.0044561 0.0104628 0.0545056 0.0076278 0.6607098 36523.62 0.2411462 NA NA NA 0.0216987 2513.664 0.2094634 1.9299725 7.892078 12.60501 NA NA NA NA NA 0.6233307 0.0685092 0.0173301 NA NA NA 0.0252765 2.3432735 0.6877957 0.0371077 0.1290979 0.1127880 0.7102951
4 401 Winston-Salem 43.22648 -0.0193413 0.2634850 0.0288748 0.0478488 0.0685591 0.4064568 30403.07 0.4511340 9.971601 0.3481184 0.5376577 0.0451547 2137.595 0.0028042 -0.0035268 7.054153 20.89541 -10.453489 NA 0.0418891 2523.8348 0.1363400 0.6357111 0.1408433 0.8333290 0.0239467 0.0580380 0.0682463 0.1903379 -0.6936674 0.3871659 0.0499824 0.2215046 0.1146979 0.5695345
5 402 Martinsville 53.63568 0.0056276 0.0069904 0.0220931 0.0410883 0.0172575 0.5862997 29943.32 0.2565523 4.503661 0.1837683 0.7311252 0.0266855 2756.848 0.3987184 -0.0275933 4.544135 18.19494 5.789057 -0.0153886 NA NA NA 0.6121494 0.1984681 0.5390288 0.0377694 0.0170314 0.0643845 0.0236150 0.3655432 0.8407629 0.0321638 0.1507503 0.1158269 0.6865350

ici nous pouvons manipuler légèrement nos données

# Example avec manipulation des données 
# Décrire ce qui se passe ici
mobility_data %>% 
  filter(!is.na(e_rank_b)) %>% 
  mutate(decile = ntile(e_rank_b, 10)) %>% 
  select(decile, e_rank_b, mig_outflow) %>% 
  group_by(decile) %>% 
  summarise(average_outflow = mean(mig_outflow)) %>% 
  kable()
decile average_outflow
1 0.0463487
2 0.0459592
3 0.0457738
4 0.0436270
5 0.0413216
6 0.0447277
7 0.0432519
8 0.0408539
9 0.0399971
10 NA

4.2 Distributions en R

L’un des exercices les plus importants lors de la modélisation statistique consiste à vérifier les hypothèses du modèle. La logique est toujours la même.

q-distribution (quantiles, ex: qnorm, qpois, qgamma etc.) Quel est le 99ème quantile d’une gamma avec shape=2 et rate=1?

qgamma(0.99, shape=2,rate=1)
## [1] 6.638352

d-distribution (densité, ex: dnorm, dpois) Tracez visuellement le 99ème quantile d’une gamma avec shape=2 et rate=1?

gamma_dist <- dgamma(seq(0,8,0.1), shape=2,rate=1)

plot(x=seq(0,8,0.1), 
     y=gamma_dist,
     type='l')
abline(v=qgamma(0.99, shape=2,rate=1), 
       col='red', lty=2)

p-distribution (densité cumulative, ex: pnorm, ppois)

Supposons que \(X \sim Gamma(2,1)\). Quel pourcentage de la masse totale se trouve à droite de 5 ?

1 - pgamma(5, shape=2, rate=1)
## [1] 0.04042768

r-distribution (échantillon, ex: rnorm, rpois)

Simulez 100 observations à partir d’un gamma avec shape=2, rate=1 et comparez-les à la distribution théorique.

set.seed(42)
sims_gamma <- rgamma(100, shape=2, rate=1)

plot(density(sims_gamma), ylim=c(0,0.38), 
     xlim=c(-1, 7), 
     main='Comparaison empirique vs. theoretique')
lines(x=seq(0,8,0.1), 
      y=gamma_dist, col='red')

4.2.1 Q-Q Plot

Pour comparer un ajustement avec une distribution théorique, on peut utiliser un diagramme quantile-quantile (qq-plot). Par exemple, on suppose souvent que le log revenue est distribué selon une normale. Comparons l’ajustement théorique et l’ajustement réel.

# A la main
# Sans transformation
inc <- (as.numeric(mobility_data$hhinc00))
inc <- (inc - mean(inc)) / sd(inc)

quantiles_norm <- qnorm(seq(0,1,0.001))
quantiles_empirical <- quantile(inc, seq(0,1,0.001))

plot(quantiles_norm,
     quantiles_empirical, 
     main='Normal Q-Q Plot', 
     ylab='Sample Quantiles',  
     xlab='Theoretical Quantiles')
abline(a = 0, b = 1, col = "red")

# Avec transformation
# A la main
# Sans transformation
log_inc <- log(as.numeric(mobility_data$hhinc00))
log_inc <- (log_inc - mean(log_inc)) / sd(log_inc)

quantiles_norm <- qnorm(seq(0,1,0.001))
quantiles_empirical <- quantile(log_inc, seq(0,1,0.001))

plot(quantiles_norm,
     quantiles_empirical, 
     main='Normal Q-Q Plot (Log transformed)', 
     ylab='Sample Quantiles',  
     xlab='Theoretical Quantiles')
abline(a = 0, b = 1, col = "red")

qqnorm(log(as.numeric(mobility_data$hhinc00)))
qqline(log(as.numeric(mobility_data$hhinc00)),
       col = 2,lwd=2,lty=2)

4.3 Exercises

  1. Supposons que \(X \sim \mathcal{N}(0,10)\). Echantillonnez 10 observations et calculez la moyenne, que nous appelons Y. Répétez l’opération ci-dessus 1000 fois et calculez la moyenne et l’écart-type de Y. A quoi vous attendez-vous et pourquoi, qu’observez-vous ?

  2. Créez un tableau ou qui montre le nombre d’observations manquantes pour chaque variable (colonne dans l’ensemble de données).

  3. Analysez les variables tuition et num_inst_pc et discutez de ce que vous voyez. Quelles sont vos conclusions ? S’agit-il d’un problème?

  4. Créez un tableau montrant le coefficient de corrélation entre la mobilité sociale e_rank_b et toutes les autres variables. Évaluez vos résultats.

  5. Séparez la variable en “faible”, “moyen” et “élevé”. Justifiez votre choix. Comparez chaque niveau avec les variables gini et cs_born_foreign. Interprétez et créez un tableau.

5 Example guidée

Dans cet exemple guidé, nous allons à nouveau travailler avec les données de Raj Chetty, Nathaniel Hendren, Patrick Kline et Emmanuel Saez.

Notre objectif est de mieux comprendre la mobilité intragénérationnelle. On dit que la mobilité intergénérationnelle est élevée si le revenu des parents (au cours de leur vie) n’est pas très lié au revenu de leurs enfants. Pour les exercices, vous pouvez utiliser n’importe quel paquet que vous jugez nécessaire mais vous devez justifier vos résultats.

Pour cela, résolvez les exercices suivants :

  1. Chargez les données, regardez les types de données, est-ce que cela a du sens ?

  2. La variable d’intérêt est e_rank_b, renommez cette colonne. Que se passe-t-il si vous utilisez un espace ” ” dans le nom ? Est-ce une bonne idée ?

  3. Visualiser la distribution de la mobilité intergénérationnelle. Donnez à vos graphiques un titre approprié et des étiquettes d’axe correctes.

  4. Examinez la relation entre la mobilité intergénérationnelle et la taille du secteur manufacturier. Quelle est votre conclusion ?

  5. Découpez la variable de la taille du secteur manufacturier en trois composantes (low-mid-high). Calculez les moyennes par groupe pour trois autres variables que vous jugez intéressantes.

  6. Ajustez une régression qui décrit l’équation suivante \(\mathbb{E}[y| \text{manufacturing}]=\hat{\beta}_0 + \hat{\beta}_1*\text{manuf_1} + \hat{\beta}_2*\text{manuf_2}\). Pourquoi exclure la dernière composante?

L’existence d’une classe moyenne importante est souvent considérée comme une bonne chose. Voir par exemple cet article, analysons cela plus en détail.

  1. Réalisez un nuage de points (“scatterplot”) avec la mobilité intergénérationnelle sur l’axe des y et la fraction de la classe moyenne sur l’axe des x. Quelles sont vos conclusions ?

  2. Exécutez un modèle linéaire qui explique la mobilité intergénérationnelle à l’aide de la fraction de la classe moyenne.

  3. Nous voulons rendre le modèle légèrement plus flexible. Effectuez 10 régressions différentes qui correspondent à des polynômes de plus en plus complexes de la fraction de la classe moyenne.

  4. Défi : Visualisez les différents résultats. Qu’arrive-t-il à la valeur \(R^2\)? Est-ce une bonne chose ?

5.1 Charger et preparer les données

# Charger l'essentiel
library(tidyverse)
library(knitr)
library(gridExtra)
library(haven) 
library(readxl)

mobility_data <- read_dta("data/mobility_across_cz.dta")
codebook <- read_excel("data/codebook.xls")

mobility_data %>% 
  head() %>% 
  kable()
index cz czname e_rank_b cs_race_bla cs_race_theil_2000 cs00_seg_inc cs00_seg_inc_pov25 cs00_seg_inc_aff75 frac_traveltime_lt15 hhinc00 gini inc_share_1perc gini99 frac_middleclass taxrate subcty_total_expenditure_pc tax_st_diff_top20 eitc_exposure ccd_exp_tot ccd_pup_tch_ratio score_r dropout_r num_inst_pc tuition gradrate_r cs_labforce cs_elf_ind_man d_tradeusch_pw_1990 frac_worked1416 mig_inflow mig_outflow cs_born_foreign scap_ski90pcm rel_tot crime_violent cs_fam_wkidsinglemom cs_divorced cs_married
0 100 Johnson City 43.58841 0.0693926 0.0576880 0.0408645 -0.0171138 0.0266398 0.3844612 38209.54 0.4193436 5.695437 0.2528637 0.5840611 0.0158913 1404.822 1.0794937 0.0000780 6.577284 17.96812 NA NA 0.0554321 1032.0113 -0.2210098 0.6529477 0.0600338 NA 0.0148774 0.0417411 0.0559070 0.1604878 NA 0.5358987 0.0482720 0.2421846 0.1392865 0.5491035
1 200 Morristown 39.64703 0.4453178 0.0749721 -0.0258754 0.0352996 0.0089049 0.3543061 25123.34 0.4867183 11.816536 0.3926104 0.3559214 0.0257373 1568.710 0.0198662 -0.0128659 4.405095 17.79943 -16.015344 0.0561887 0.0379247 970.0455 -0.0416078 0.5273348 0.2086723 0.7670213 0.0223866 0.0181590 0.0137850 0.0232165 -0.6127989 0.6469429 0.0316044 0.3499764 0.1039057 0.5426201
2 301 Middlesborough 46.93892 0.0151223 0.1100580 0.0954185 0.0139200 0.0527287 0.4793166 32701.43 0.3071970 11.810008 0.2077348 0.6675747 0.0116976 3115.946 -0.0046168 0.0343010 6.609784 17.85102 -9.137633 0.0246227 0.0634999 1395.0190 -0.1716805 0.7066224 0.0362908 0.0274560 0.0274638 0.0309692 0.0421808 0.1318125 -2.0635338 0.4222386 0.0511852 0.1831544 0.1402021 0.6224286
3 302 Knoxville NA -0.0004112 -0.0044561 0.0104628 0.0545056 0.0076278 0.6607098 36523.62 0.2411462 NA NA NA 0.0216987 2513.664 0.2094634 1.9299725 7.892078 12.60501 NA NA NA NA NA 0.6233307 0.0685092 0.0173301 NA NA NA 0.0252765 2.3432735 0.6877957 0.0371077 0.1290979 0.1127880 0.7102951
4 401 Winston-Salem 43.22648 -0.0193413 0.2634850 0.0288748 0.0478488 0.0685591 0.4064568 30403.07 0.4511340 9.971601 0.3481184 0.5376577 0.0451547 2137.595 0.0028042 -0.0035268 7.054153 20.89541 -10.453489 NA 0.0418891 2523.8348 0.1363400 0.6357111 0.1408433 0.8333290 0.0239467 0.0580380 0.0682463 0.1903379 -0.6936674 0.3871659 0.0499824 0.2215046 0.1146979 0.5695345
5 402 Martinsville 53.63568 0.0056276 0.0069904 0.0220931 0.0410883 0.0172575 0.5862997 29943.32 0.2565523 4.503661 0.1837683 0.7311252 0.0266855 2756.848 0.3987184 -0.0275933 4.544135 18.19494 5.789057 -0.0153886 NA NA NA 0.6121494 0.1984681 0.5390288 0.0377694 0.0170314 0.0643845 0.0236150 0.3655432 0.8407629 0.0321638 0.1507503 0.1158269 0.6865350
str(mobility_data)
## tibble [495 × 39] (S3: tbl_df/tbl/data.frame)
##  $ index                      : num [1:495] 0 1 2 3 4 5 6 7 8 9 ...
##   ..- attr(*, "format.stata")= chr "%12.0g"
##  $ cz                         : num [1:495] 100 200 301 302 401 402 500 601 602 700 ...
##   ..- attr(*, "format.stata")= chr "%12.0g"
##  $ czname                     : chr [1:495] "Johnson City" "Morristown" "Middlesborough" "Knoxville" ...
##   ..- attr(*, "format.stata")= chr "%22s"
##  $ e_rank_b                   : num [1:495] 43.6 39.6 46.9 NA 43.2 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_race_bla                : num [1:495] 0.069393 0.445318 0.015122 -0.000411 -0.019341 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_race_theil_2000         : num [1:495] 0.05769 0.07497 0.11006 -0.00446 0.26349 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc               : num [1:495] 0.0409 -0.0259 0.0954 0.0105 0.0289 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc_pov25         : num [1:495] -0.0171 0.0353 0.0139 0.0545 0.0478 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc_aff75         : num [1:495] 0.02664 0.0089 0.05273 0.00763 0.06856 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_traveltime_lt15       : num [1:495] 0.384 0.354 0.479 0.661 0.406 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ hhinc00                    : num [1:495] 38210 25123 32701 36524 30403 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gini                       : num [1:495] 0.419 0.487 0.307 0.241 0.451 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ inc_share_1perc            : num [1:495] 5.7 11.82 11.81 NA 9.97 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gini99                     : num [1:495] 0.253 0.393 0.208 NA 0.348 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_middleclass           : num [1:495] 0.584 0.356 0.668 NA 0.538 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ taxrate                    : num [1:495] 0.0159 0.0257 0.0117 0.0217 0.0452 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ subcty_total_expenditure_pc: num [1:495] 1405 1569 3116 2514 2138 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ tax_st_diff_top20          : num [1:495] 1.07949 0.01987 -0.00462 0.20946 0.0028 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ eitc_exposure              : num [1:495] 0.000078 -0.012866 0.034301 1.929972 -0.003527 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ ccd_exp_tot                : num [1:495] 6.58 4.41 6.61 7.89 7.05 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ ccd_pup_tch_ratio          : num [1:495] 18 17.8 17.9 12.6 20.9 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ score_r                    : num [1:495] NA -16.02 -9.14 NA -10.45 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ dropout_r                  : num [1:495] NA 0.0562 0.0246 NA NA ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ num_inst_pc                : num [1:495] 0.0554 0.0379 0.0635 NA 0.0419 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ tuition                    : num [1:495] 1032 970 1395 NA 2524 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gradrate_r                 : num [1:495] -0.221 -0.0416 -0.1717 NA 0.1363 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_labforce                : num [1:495] 0.653 0.527 0.707 0.623 0.636 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_elf_ind_man             : num [1:495] 0.06 0.2087 0.0363 0.0685 0.1408 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ d_tradeusch_pw_1990        : num [1:495] NA 0.767 0.0275 0.0173 0.8333 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_worked1416            : num [1:495] 0.0149 0.0224 0.0275 NA 0.0239 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ mig_inflow                 : num [1:495] 0.0417 0.0182 0.031 NA 0.058 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ mig_outflow                : num [1:495] 0.0559 0.0138 0.0422 NA 0.0682 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_born_foreign            : num [1:495] 0.1605 0.0232 0.1318 0.0253 0.1903 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ scap_ski90pcm              : num [1:495] NA -0.613 -2.064 2.343 -0.694 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ rel_tot                    : num [1:495] 0.536 0.647 0.422 0.688 0.387 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ crime_violent              : num [1:495] 0.0483 0.0316 0.0512 0.0371 0.05 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_fam_wkidsinglemom       : num [1:495] 0.242 0.35 0.183 0.129 0.222 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_divorced                : num [1:495] 0.139 0.104 0.14 0.113 0.115 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_married                 : num [1:495] 0.549 0.543 0.622 0.71 0.57 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"

Le format est celui de stata, mais la plupart des variables sont numériques, cela semble correct.

# Verifier données manquantes
for (column_ in names(mobility_data)){
  print(
    paste(
      column_, 
      sum(is.na(mobility_data[, column_])), 
      sep=': '
    )
  )
}
## [1] "index: 0"
## [1] "cz: 0"
## [1] "czname: 0"
## [1] "e_rank_b: 19"
## [1] "cs_race_bla: 0"
## [1] "cs_race_theil_2000: 0"
## [1] "cs00_seg_inc: 0"
## [1] "cs00_seg_inc_pov25: 0"
## [1] "cs00_seg_inc_aff75: 0"
## [1] "frac_traveltime_lt15: 0"
## [1] "hhinc00: 0"
## [1] "gini: 0"
## [1] "inc_share_1perc: 19"
## [1] "gini99: 19"
## [1] "frac_middleclass: 19"
## [1] "taxrate: 1"
## [1] "subcty_total_expenditure_pc: 2"
## [1] "tax_st_diff_top20: 0"
## [1] "eitc_exposure: 0"
## [1] "ccd_exp_tot: 7"
## [1] "ccd_pup_tch_ratio: 20"
## [1] "score_r: 21"
## [1] "dropout_r: 89"
## [1] "num_inst_pc: 100"
## [1] "tuition: 104"
## [1] "gradrate_r: 107"
## [1] "cs_labforce: 0"
## [1] "cs_elf_ind_man: 0"
## [1] "d_tradeusch_pw_1990: 14"
## [1] "frac_worked1416: 19"
## [1] "mig_inflow: 10"
## [1] "mig_outflow: 10"
## [1] "cs_born_foreign: 0"
## [1] "scap_ski90pcm: 14"
## [1] "rel_tot: 0"
## [1] "crime_violent: 17"
## [1] "cs_fam_wkidsinglemom: 0"
## [1] "cs_divorced: 0"
## [1] "cs_married: 0"

Le nombre de points de données manquants ne semble pas non plus trop inquiétant. Il pourrait y avoir une certaine corrélation entre les points manquants mais nous pouvons ignorer cela pour le moment.

5.2 Renommer une variable

# renommer avec espace
print(names(mobility_data))
##  [1] "index"                       "cz"                         
##  [3] "czname"                      "e_rank_b"                   
##  [5] "cs_race_bla"                 "cs_race_theil_2000"         
##  [7] "cs00_seg_inc"                "cs00_seg_inc_pov25"         
##  [9] "cs00_seg_inc_aff75"          "frac_traveltime_lt15"       
## [11] "hhinc00"                     "gini"                       
## [13] "inc_share_1perc"             "gini99"                     
## [15] "frac_middleclass"            "taxrate"                    
## [17] "subcty_total_expenditure_pc" "tax_st_diff_top20"          
## [19] "eitc_exposure"               "ccd_exp_tot"                
## [21] "ccd_pup_tch_ratio"           "score_r"                    
## [23] "dropout_r"                   "num_inst_pc"                
## [25] "tuition"                     "gradrate_r"                 
## [27] "cs_labforce"                 "cs_elf_ind_man"             
## [29] "d_tradeusch_pw_1990"         "frac_worked1416"            
## [31] "mig_inflow"                  "mig_outflow"                
## [33] "cs_born_foreign"             "scap_ski90pcm"              
## [35] "rel_tot"                     "crime_violent"              
## [37] "cs_fam_wkidsinglemom"        "cs_divorced"                
## [39] "cs_married"
names(mobility_data)[4] <- "mobilite intergen"

# mobility_data$`mobilite intergen` -> On a toujours besoin des symboles `. Cela semble encombrant. 
# Essayez d'utiliser des tirets bas (_) ou des points (.)

names(mobility_data)[4] <- "mobilite_intergen"
# mobility_data$mobilite_intergen -> C'est mieux

5.3 Analyse visuelle de la mobilité

Nous allons analyser la densité en utilisant le lissage par noyau, un histogramme et la fonction de distribution cumulative.

hist(mobility_data$mobilite_intergen, 
     main='Intergenerational Mobility Scores', 
     xlab='Value', 
     ylab='Density', 
     probability = T, 
     breaks = seq(27, 67, 3))

lines(density(mobility_data$mobilite_intergen[!is.na(mobility_data$mobilite_intergen)]), 
      col='red', 
      lwd=2, 
      lty=2)

legend(55,0.07,
       legend=c("Density"),
       col=c("red"),
       lty=2,
       cex=0.8)

quantiles_ <- quantile(mobility_data$mobilite_intergen, 
                       probs = seq(0,1,0.01), 
                       na.rm = T)

plot(x=quantiles_, 
     y=seq(0,1,0.01), 
     type='l', 
     lwd=2, 
     main='ECDF', 
     ylab='', 
     xlab='Value for Mobility')

5.4 Analyse visuelle: Manufacturing vs Mobilité

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

# Model simple pour la relation lineaire
summary(lm(mobilite_intergen ~ cs_elf_ind_man, data=mobility_data))
## 
## Call:
## lm(formula = mobilite_intergen ~ cs_elf_ind_man, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5226  -3.6620  -0.8982   3.1683  17.9719 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     46.9024     0.5867  79.946  < 2e-16 ***
## cs_elf_ind_man -17.7064     3.1306  -5.656 2.68e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.482 on 474 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.06322,    Adjusted R-squared:  0.06125 
## F-statistic: 31.99 on 1 and 474 DF,  p-value: 2.683e-08
plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

abline(a=46.9024,
       b=-17.7064,
       col='red', 
       lwd=2, 
       lty=2)

legend(0.35,60,
       legend=c("Linear Fit"),
       col=c("red"),
       lty=2,
       cex=0.8)

5.5 Analyse descriptif (secteur manufacturier)

Nous pouvons trouver de bons seuils en regardant le graphique ci-dessus ou en utilisant les quantiles. Voici l’approche “visuelle”. Lorsque nous coupons une variable, nous pouvons introduire des non-linéarités dans la relation (voir l’exercice suivant)

mobility_data$manuf_cut <- cut(mobility_data$cs_elf_ind_man, 
                               breaks = c(0,0.1,0.3, 1), 
                               labels = c('low', 'mid', 'high'))

mobility_data %>% 
  group_by(manuf_cut) %>% 
  summarise(
    mean_labforce = mean(cs_labforce, na.rm=T), 
    mean_growth_imports = mean(d_tradeusch_pw_1990, na.rm=T), 
    mean_mig_inflow = mean(mig_outflow, na.rm=T)
  )
## # A tibble: 3 × 4
##   manuf_cut mean_labforce mean_growth_imports mean_mig_inflow
##   <fct>             <dbl>               <dbl>           <dbl>
## 1 low               0.633               0.289          0.0458
## 2 mid               0.641               1.30           0.0427
## 3 high              0.667               2.73           0.0399

5.6 Model simple

model_lineaire_simple = lm(mobilite_intergen ~ manuf_cut, 
                           data=mobility_data)

summary(model_lineaire_simple)
## 
## Call:
## lm(formula = mobilite_intergen ~ manuf_cut, data = mobility_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.271  -3.847  -0.658   2.952  17.605 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    45.9575     0.5388  85.293  < 2e-16 ***
## manuf_cutmid   -2.3677     0.6141  -3.856 0.000131 ***
## manuf_cuthigh  -6.3979     1.2444  -5.142    4e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.495 on 473 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.06088,    Adjusted R-squared:  0.05691 
## F-statistic: 15.33 on 2 and 473 DF,  p-value: 3.534e-07
plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

x <- seq(0, 0.45, 0.01)
fx <- 45.9575 + (x > 0 & x <=0.1) * 0 +
      (x > 0.1 & x <= 0.3) * -2.3677 +
      (x > 0.3 & x <= 2.119) *  -6.3979
lines(x, fx, 
      col='red', 
      lty=2, 
      lwd=2)

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_simple <- predict(model_lineaire_simple, 
                              mobility_data)

points(mobility_data$cs_elf_ind_man, 
       predictions_simple, 
       col='red')

model_lineaire_interaction = lm(mobilite_intergen ~ manuf_cut + cs_elf_ind_man + manuf_cut*cs_elf_ind_man, 
                           data=mobility_data)

summary(model_lineaire_interaction)
## 
## Call:
## lm(formula = mobilite_intergen ~ manuf_cut + cs_elf_ind_man + 
##     manuf_cut * cs_elf_ind_man, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.9773  -3.6249  -0.6577   2.7890  17.3450 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    42.921      1.983  21.640   <2e-16 ***
## manuf_cutmid                    3.175      2.252   1.410   0.1592    
## manuf_cuthigh                  -5.908     10.728  -0.551   0.5821    
## cs_elf_ind_man                 44.485     27.973   1.590   0.1124    
## manuf_cutmid:cs_elf_ind_man   -57.912     28.507  -2.031   0.0428 *  
## manuf_cuthigh:cs_elf_ind_man  -37.336     40.598  -0.920   0.3582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.463 on 470 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.07768,    Adjusted R-squared:  0.06787 
## F-statistic: 7.917 on 5 and 470 DF,  p-value: 3.562e-07
plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_interact <- predict(model_lineaire_interaction, 
                              mobility_data)

points(mobility_data$cs_elf_ind_man, 
       predictions_interact, 
       col='red')

5.7 Model “Classe moyenne”

plot(x=mobility_data$frac_middleclass, 
     y=mobility_data$mobilite_intergen, 
     main='Middle Class and Intergen. Mobility', 
     ylab='Intergen. Mobility', 
     xlab='Fraction Middle class')

model_middle_class <- lm(mobilite_intergen ~ frac_middleclass, 
                         data=mobility_data)

summary(model_middle_class)
## 
## Call:
## lm(formula = mobilite_intergen ~ frac_middleclass, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.9940  -2.9766  -0.5619   2.7327  15.9744 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        18.986      1.381   13.74   <2e-16 ***
## frac_middleclass   45.411      2.491   18.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.343 on 474 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.4121, Adjusted R-squared:  0.4109 
## F-statistic: 332.3 on 1 and 474 DF,  p-value: < 2.2e-16
plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_middleclass <- predict(model_middle_class, 
                              mobility_data)

points(mobility_data$frac_middleclass,
       predictions_middleclass, 
       col='red')

# Analyse des residues

residuals_ <- mobility_data$mobilite_intergen - predictions_middleclass
plot(mobility_data$frac_middleclass, 
     residuals_, 
     ylab='Residuals', 
     xlab='Middle class')

mobility_data_nomissings <- mobility_data[!is.na(mobility_data$frac_middleclass), ]

model_middle_class_poly <- lm(mobilite_intergen ~ poly(frac_middleclass, 2), 
                              data=mobility_data_nomissings)

summary(model_middle_class_poly)
## 
## Call:
## lm(formula = mobilite_intergen ~ poly(frac_middleclass, 2), data = mobility_data_nomissings)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0057  -2.9207  -0.3494   2.4430  15.5096 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 43.9039     0.1881 233.464  < 2e-16 ***
## poly(frac_middleclass, 2)1  79.1684     4.1029  19.296  < 2e-16 ***
## poly(frac_middleclass, 2)2  31.2715     4.1029   7.622 1.38e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.103 on 473 degrees of freedom
## Multiple R-squared:  0.4764, Adjusted R-squared:  0.4742 
## F-statistic: 215.2 on 2 and 473 DF,  p-value: < 2.2e-16
plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_middleclass_poly <- predict(model_middle_class_poly,
                                        mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ])

lines(mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ]$frac_middleclass,
       predictions_middleclass_poly, 
       col='red')

5.8 Challenge

# List pour les resultats
results_list <- list()

# Loop
for (length_ in seq(1,10,1)){
  reg_ <- lm(mobilite_intergen ~ poly(frac_middleclass, length_), 
             data=mobility_data_nomissings)
  
  print(paste('R-Squared: ', summary(reg_)$r.squared))
  
  predictions_ <- predict(reg_,
                          mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ])
  
  results_list[[length_]] <- predictions_

}
## [1] "R-Squared:  0.412134418340672"
## [1] "R-Squared:  0.476437415437812"
## [1] "R-Squared:  0.477700187650047"
## [1] "R-Squared:  0.477727823275877"
## [1] "R-Squared:  0.478923076244006"
## [1] "R-Squared:  0.48017980394324"
## [1] "R-Squared:  0.48426795442976"
## [1] "R-Squared:  0.487546350228303"
## [1] "R-Squared:  0.487606950061351"
## [1] "R-Squared:  0.489193655479349"
plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

for(i in seq(1,10,1)){
  
  predictions_middleclass_poly <- results_list[[i]]
  
  lines(mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ]$frac_middleclass,
         predictions_middleclass_poly, 
         col='red')
  
}